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Self-propelled rods exhibit a novel phase-separated state characterized by the 
presence of active stresses and the ejection of polar clusters 
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We study collections of self-propelled rods (SPR) moving in two dimensions for packing fractions 
less than or equal to 0.3. We find that in the thermodynamical limit the SPR undergo a phase 
transition between a disordered gas and a novel phase-separated system state. Interestingly, (global) 
orientational order patterns - contrary to what has been suggested - vanish in this limit. In the 
found novel state, the SPR self-organize into a highly dynamical, high-density, compact region - 
which we call aggregate - which is surrounded by a disordered gas. Active stresses build inside 
aggregates as result of the combined effect of local orientational order and active forces. This leads 
to the most distinctive feature of these aggregates: constant ejection of polar clusters of SPR. This 
novel phase-separated state represents a novel state of matter characterized by large fluctuations 
in volume and shape, related to mass ejection, and exhibits positional as well as orientational local 
order. SPR systems display new physics unseen in other active matter systems. 

PACS numbers: 05.65.+b, 87.18.Gh, 87.18.Ed,47.54.-r 


I. INTRODUCTION 

Self-organized patterns of self-propelled entities - from 
animals to synthetic active particles - are often believed 
to be the result of a velocity alignment mechanism that 
regulates the interaction among the moving entities [l|- 
0. Such a velocity alignment mechanism can result, for 
instance, from hydrodynamic interactions 0^, electro¬ 
static forces moving molecular motors in arrays of 
microtubules [8|, loj , or from inelastic bouncing in driven 
granular particles fiol-ll^. There is, in addition, another 
simple and rather general way of producing a (velocity) 
alignment mechanism [T^: the combined effect of steric 
interactions among elongated objects and self-propulsion 
in a dissipative medium, which has been shown to lead 
to interesting collective phenomena [iMl. This mech¬ 
anism is at work in a broad range of active systems: 
gliding bacteria [TgI, [17|. driven granular rods [TsI, [T^ . 
chemically-driven rods 23 , [2l| , and it has been recently 
argued that also - neglecting hydrodynamic effects over 
steric effects - in swimming bacteria |22l-[^ and motility 
assays [25|, [26|. 

Here, we focus on the large-scale physical properties of 
collections of self-propelled rods (SPR) for packing frac¬ 
tions T] less than or equal to 0.3. We find that the com¬ 
bined effect of steric repulsive forces and active forces lead 
to a complex interplay between local orientational order 
and active stresses. As result of such interplay, SPR ex¬ 
hibit, for large enough system sizes, new physics unseen 
in other active matter systems. The novel phenomena 
reported here passed unnoticed in previous SPR stud¬ 
ies that were performed with either small system sizes 
or lack a finite size study pRijlsI. 1^ . In par¬ 

ticular, we provide strong evidence that global orienta¬ 
tional order patterns (see Fig. [1]), suggested to exist for 


1 polar orde r ^ nematic order ^ average cluster size 

rri=o.i5 n 

0.8 ti = 0-3 - 0.8 0.8 • 



2 4 ^6 8 10 2 4 ^6 8 10 2 4 ^6 8 10 

l\ i\ K 



FIG. I: Transition from disorder to order and phase separa¬ 
tion at small system sizes. Top, from left to right: polar order 
parameter Si, nematic order parameter S' 2 , and average clus¬ 
ter size {m) over system size N, respectively, as function of 
the particle aspect ratio k for various packing fractions 77 . 
N = 10000. Bottom: A simulation snapshot for = 10 and 
T] = 0.3, with Si = 0.5, S 2 = 0.26, {m)/N = 0.14. Notice 
that clusters are polar. The orientation of rods is color coded 
as indicated in the bottom left panel. This convention is also 
used in Figs. [Hand [5l 


T] < 0.3 ( 23 , [ 33 I, vanish in the thermodynamical limit. 
More importantly, we find that SPR undergo a phase 
transition between a disordered gas and a novel phase- 
separated system state characterized by the presence of 
active stresses and the ejection of polar clusters of SPR. 

The paper is organized as follows. We start out by 
introducing the self-propelled rod model in Sec. [Ill The 
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statistical features of the phase-transition in finite sys¬ 
tems is studied in Sec. mil A finite-size analysis is pre¬ 
sented in Sec. HYi The remarkable dynamics displayed 
by aggregates in large systems is analyzed in Sec. [Vl In 
Sec.ini we discuss the obtained results. 


II. MODEL DEFINITION 


Our model consists of N SPR moving in a two- 
dimensional space of linear size L with periodic boundary 
conditions. Each rod is driven by an active stress/force 
F that is applied along the long axis of the particle. In¬ 
teractions among rods are modeled through a repulsive 
potential, which we denote, for the Tth particle, by Ui. 
The substrate where the rods move acts as a momentum 
sink. There are three friction drag coefficients, C||, 
and ( 0 , which correspond to the drags experienced by the 
rods as the rod moves along the long axis, perpendicular 
to it, or as it rotates, respectively. In the over-damped 
limit, the equations of motion of the Tth rod are given, 
as in [Il,by: 


Xi = ^l[-VUi + FV{0i) + (Ti{t)] 


Ce [ 


dUi ,,; 


( 1 ) 

( 2 ) 


where the dot denotes a temporal derivative, corre¬ 
sponds to the position of the center of mass and Oi the 
orientation of the long axis of the rod. In Eq. o, M 
is the mobility tensor defined as fj. = Cy ^V(0j)V(0j) + 

Cl^V_L(6»i)V_L(6»i), with V(6») = {cos(6»), sin(6»)) and 
V_l( 6>) such that V(6 >).V_l( 6>) = 0. We use C|| = 10, 
= 25, Co = 2, and F = 0.4, which corresponds to 
an active speed vq = F/C\\ = 0.04. Other friction coef¬ 
ficient values, as well as drag friction models, have been 
also tested as detailed in [^, obtaining the same qual¬ 
itative results. Eor details about how to compute drag 
friction coefficients we refer the reader to (sH, [39| . The 
temporal evolution of the orientation of the rod, given 
by Eq. ([2]), results from the torque — generated by 
the interactions and no active torque is present. Eqs. 
o and m are subject to fluctuations through the terms 
cri{t) and which correspond to delta-correlated vec¬ 
torial and scalar noise, respectively. If these fluctuations 
are of thermal origin, it can be shown that the “pas¬ 
sive” diffusion coefficient Dp resulting from cri (i.e. the 
diffusion for F = 0) is negligible compared to the ac¬ 
tive diffusion coefficient given by Da (corresponding to 
F > 0 and ai = 0); for details see Appendix [Al Since 
Dp/Da ^ 1, for simplicity we neglect cri and specify in 
Eq. (|2]) {Ci{t)) = 0 and {Ci{t)Cj{t')) = 2De6ij6{t - t'), 
with Dq = 2.52 X 10“^. The interactions among the rods 
are modeled by a soft-core potential that penalizes par¬ 
ticle overlapping. Eor the Fth rod, the potential takes 
the form: Ui = U{'Ki,0i) = ^ where Uij 

denotes the repulsive potential interaction between the 
i-th and j-th rod, both of length £ and width re, such 



FIG. 2: Interactions in the SPR model. Top: Sketch of two 
interacting rods in the SPR model. Forces and torques re¬ 
sulting from the interaction are shown in the top-right panel. 
See text for explanations. Bottom) Chronological snapshots 
of a collision between two rods. Notice that even though the 
interaction is exclusively repulsive, it leads to an effective ve¬ 
locity alignment (an effective attraction). A movie of this 
interaction is provided in a. 


that Uij = u{'Ki — ^j^Oi — Oj). The rods are represented 
as straight chains of n disks of diameter re, as imple¬ 
mented in whose centers are separated at dis¬ 

tance A = rc/S. We notice that results obtained with 
SPR represented by disk-chains are qualitatively identi¬ 
cal to those produced with the original SPR model in¬ 
troduced in [1^ if and only if A C 2n;; for more details 
see [i^. Using this implementation, Uij can be expressed 
as Ui^j = ^ where u/ f is the potential between 

disk a of the Fth rod and disk P of the j-th rod, which 
here we assume to be given by a harmonic repulsive po¬ 
tential: 14= Co 1 ^ — : for < re and zero 

otherwise, where is the distance between the centers 
of the disks and Co = 200. 

Fig. [2] illustrates the implementation of the interac¬ 
tion between two rods. Rod i is propelled forward by 
its active force F\{0i)^ while it is pushed away by rod 
j through the interaction force —VUi. The interaction 
with rod j also leads to a torque, given by — This 
torque, together with the over-damped dynamics, lead 
to an effective alignment of the velocity of the rods as 
shown in the bottom panels of Fig. [2] and previously de¬ 
scribed in M- The snapshots show, in chronological 
order, a collision event between two rods: The rods are 
moving in different directions (indicated by the arrows) 
before the collision (first snapshot). They start colliding 
(second snapshot). The steric forces and torques result¬ 
ing from the interaction lead to an effective alignment 
(third snapshot). As result, both rods end up moving in 
roughly the same direction and stay close to each other 
(last snapshot) - without requiring attractive force. No¬ 
tice that interestingly the described alignment process is 
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similar to the one predicted to occur using kinetic the¬ 
ory in model for ordering of microtubules mediated by 
molecular motors in a planar geometry [I, ^ . 

III. TRANSITION FROM DISORDER TO 
ORDER AND PHASE SEPARATION IN SMALL 
FINITE SYSTEMS 

We explore the parameter space of the model us¬ 
ing parallelized numerical simulations implementing an 
Euler-Maruyama scheme. As control parameter, we use 
the rod aspect ratio k = ijw and keep constant the par¬ 
ticle area a = i x w = 0.1 diS well as all the other pa¬ 
rameters. We perform the analysis for several values of 
77. Notice that 77 can be also used as control parameter. 
The macroscopic patterns are characterized by their level 
of orientational order through 

Sq = {Sg{t))t = (|(exp(^g6»^(^)))i|)^, (3) 

where the averages run over the number of particles and 
time. Polar order corresponds to g = 1 and nematic or¬ 
der to g = 2 . Phase separation is monitored by looking 
at the ratio between average cluster size (m) and sys¬ 
tem size A, where (m) = with p{m) the 

(weighted steady state) cluster size distribution. p{m) is 
defined as the time average (after an initial transient) of 
the instantaneous cluster size distribution 

m nm{t) 

p{m,t) = - - -, (4) 

where nm{t) is the number of clusters of mass m present 
in the system at time t. Notice that the normalization 
of this distribution is ensured since N = m 

Clusters are collections of interconnected rods, where any 
two rods are considered as connected if they are separated 
by a distance equal to or smaller than 2w^ implemented 
as the minimum distance between the centers of the disks 
that form each rod. 

Varying the aspect ratio k, while keeping fixed the 
packing fraction 77 = aNjL?^ we observe a clear phase 
transition with ^Si, and {m)lN taking off above a 
critical /^c value as shown in Fig. [TJ We notice that Kc 
decreases when 77 is increased. Below /^c, we observe a gas 
phase characterized by the absence of orientational or¬ 
der and an exponential cluster size distribution such that 
{m)/N ^ 0{N~^). For k > Kc the system undergoes a 
symmetry breaking as observed previously in [l5|, [^, . 

The emerging order is polar as evidenced by the behavior 
of Si. We recall that in the presence of polar order, S 2 
is slaved to Si. The behavior of {m)/N, right panel in 
Fig. [U indicates that the system starts to spontaneously 
self-segregate for n > tic- Here, we find that the onsets 
of orientational order and phase separation coincide and 
share the same critical point, as predicted using a simple 
kinetic model for the clustering process 0 ; due to the 
effective velocity alignment large polar clusters emerge. 



FIG. 3: Finite size scaling of the polar order parameter Si 
and average cluster size with respect to system size {m) /N for 
aspect ratio tc = 10 and several packing fractions 77 . Notice 
that, for 77 > 0.15, the system becomes disordered as the 
system size N is increased, while remaining phase-separated. 


which in turn lead to macroscopic polar order. Notice 
that in an equilibrium system of (hard) rods (i.e. F = 0) 
for 72< 0.3 and 1 < k < 10^ according to De las Her as et 
al. [^, we should observe only an homogeneous disor¬ 
dered phase for this range of parameters. This indicates 
that the observed phase transition requires F > 0, i.e. 
the active motion of the rods. 


IV. FINITE SIZE STUDY: ABSENCE OF 
GLOBAL ORDER IN THE 
THERMODYNAMICAL LIMIT 

We performed a finite size study, by increasing simul¬ 
taneously N and L while keeping the packing fraction 77 
and all other parameters constant. Fig.[3]shows the scal¬ 
ing of the (global) polar order parameter Si and average 
cluster size with respect to system size {m)/N for aspect 
ratio hz = 10 and several packing fractions 77 . At low 77 
values, i.e. for 77 < 0.075, we are in the situation k < Kc 
(we recall that the critical Kc value depends on 77 ). We 
showed in Sec. IHII that for K < hZc the system is not phase- 
separated and does not exhibit orientational order. As 
expected, the scaling of Si with N shows that Si oc 
with a = 1 / 2 , which means that the system is fully dis¬ 
ordered. In addition, we observe that oc A“^, with 
yd = 1 , which indicates that there is a well-defined charac¬ 
teristic cluster size for the system (that is independent of 
A) and consequently the system is not phase-separated. 

At large 77 values, i.e. for 77 > 0.15 and hz > tZc^ we ob¬ 
serve phase separation and (global) polar order for small 
finite systems (A = 10000), Sec. [Till The finite size study 
shows that, for k > hZc, does not decrease (asymp¬ 

totically) with A. Moreover, for large values of A, we 
even observe an increase. This indicates that {m) is at 
least proportional to A, which implies that the system is 
phase-separated in the thermodynamical limit as well as 
in finite systems. At the level of the orientational order 
parameter Fi, we observe an abrupt change in scaling of 
Si with A. For N < (e.g., ^ 20000) for 77 = 0.3 

and hz = 10, Si is high i.e. the system displays global po¬ 
lar order . On the other hand, for A > A^, while the sys¬ 
tem remains phase-separated. Si sharply decreases with 
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N. In short, the finite size study reveals that, although 
phase separation does take place in the thermodynamical 
limit as well as in finite systems, the phase transition to 
an orientationally ordered phase, described in Sec. Iml is 
observed only for small finite systems. Global order pat¬ 
terns are not present in the thermodynamical limit, with 
the phase transition occurring, in this limit, between a 
disordered gas and a phase-separated state with no global 
orientational order. 

The reason for observing non-vanishing global order 
in small systems is the presence of few giant polar clus¬ 
ters as illustrated by the simulation snapshot in Fig. [TJ 
Such giant polar clusters can become so big and elon¬ 
gated that they can can even percolate the system, as 
shown in panel i) of Fig. 01 We refer to such polar perco¬ 
lating structures as bands. Inside bands, rods are densely 
packed, point into the same direction, and exhibit posi¬ 
tional order. Notice that these bands are distinct from 
the bands observed in the Vicsek model, which are elon¬ 
gated in the direction orthogonal to the moving direc¬ 
tion of the particles [i2|. The observed polar bands are 
also fundamentally different from those observed in mod¬ 
els of idealized SPR, where the point-like self-propelled 
particles form nematic bands, inside which 50% of the 
particles move in one direction and 50% in the opposite 
one [ 29 I . More importantly, our finite size study indicates 
that the polar patterns observed in SPR are a finite size 
effect that disappear for large enough systems. In short, 
several of the phases reported for r] < 0.3 in previous SPR 
works [ 22 I, [ 33 I such as the so-called swarming phase and 
the bio-turbulence phase vanish in the thermodynamical 
limit. 

The abrupt change in scaling of Si with N in Fig. [3] 
suggests that above the crossover system size the po¬ 
lar structures are no longer stable. Arguably, the decay 
in Si with N is due to the fact that rods inside polar 
clusters are densely packed and hold fixed positions, not 
being able to exchange neighbors in contrast to other ac¬ 
tive systems [28|, [29|, 1^ . [sj, EH . In the co-moving frame 
that moves with the cluster, we have a two-dimensional 
system of particles interacting locally and subject to fluc¬ 
tuations. Assuming that in this scenario we can apply 
the Mermin-Wagner theorem jHH , long-range order is not 
possible and for sufficiently big clusters defects in the ori¬ 
entation of the rods should emerge. If such defects are 
present in a cluster/band, the velocity field of the po¬ 
lar structure will be necessarily unstable (see also Ap¬ 
pendix!^. The instability of polar structures is evident 
by looking at the behavior of bands with N, Fig. jU The 
top panel of this figure shows the finite size scaling of 
TaggjTtot^ the total time Tagg the system spends in 
the aggregate phase with respect to the total simulation 
time Ttot‘ Note that the computation of Tagg implies 
looking for all events where an aggregate emerged in the 
system, accumulating the time each aggregate lived. We 
observe that Tagg increases with N^ in such a way that 
TaggjTtot ^ las ^ 00 . This means that the proba¬ 
bility of observing the system in an aggregate phase also 



FIG. 4: The orientationally ordered phase-separated state 
observed for > /^c in small systems {N N^) becomes in¬ 
stable when N is increased. In very big systems {N ^ W) we 
observe only aggregates (which correspond to an orientation¬ 
ally disordered phase-separated state). Top left: Evolution of 
TaggiTtot with A, whcrc Tagg is the time the system spent 
in the aggregate phase and Ttot - here, Ttot = 10 ^ - is the 
total simulation time. Second line, left: the polar order Si as 
function of time for N - here, A* ~ 10^. Large values 

of Si correspond to the system being highly ordered, typically 
due to the formation of a highly ordered band (panel i), while 
low values of Si correspond to the formation of an aggregate 
(panel ii). Third line: Histograms of the polar order param¬ 
eter p{Si) for A = 5000, A = 10000 - A. and A = 20000. 
Bottom panel: Total elastic energy Utot = 
system as function of time t, for N ^ N^. High values of Utot 
correspond to the presence of an aggregate (inset on the left 
panel), while low values of Utot are related to the presence of 
polar structures such as bands. On both snapshots, the color 
of each rod indicates its interaction potential Up. blue (red) 
color indicates small (large) Ui values. The dashed black cir¬ 
cle in the snapshot corresponding to an aggregate provides 
an idea of the boundary shell of the aggregate. For movies, 
see EH- Simulations correspond to = 4 and 77 = 0.3. 


increases with A. For small system sizes A ^ A^ we ob¬ 
serve moving clusters and bands. Large polar structures 
such as bands form, remain in the system for quite some 
time, and then quickly break and reform, typically adopt¬ 
ing a new orientation. The corresponding histogram of 
global polar order - p{Si) - in Fig. [4] (A = 5000) is uni- 
modal with a peak at large values oi Si. As A ^ A*, 
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bands survive for relatively short periods of time, and 
quickly bend and break. Interestingly, at such large sys¬ 
tem sizes other macroscopic structures start to frequently 
emerge. These new macroscopic structures - which we 
refer to as aggregates - are formed by polar clusters of 
rods that exert stresses on each other and exhibit van¬ 
ishing polar order, see panel ii) of Fig. 01 In summary, 
for N ^ the system continuously transitions between 
highly ordered phases - e.g. phases with either a few 
giant polar clusters or a band - and aggregates, as il¬ 
lustrated in Fig. m The corresponding histogram of Si 
{N = 10000) is bimodal with a peak at large values of 
Si, corresponding to the polar structures, and another 
peak at very small values of Si, corresponding to aggre¬ 
gates. As the system size is increased further, i.e. for 
N ^ N^, we observe that the corresponding histogram 
of Si {N = 20000) becomes again unimodal, but the peak 
is now at very small values of Si, and corresponds to the 
presence of aggregates. In short, bands and polar phases 
disappear in the thermodynamical limit, while the ag¬ 
gregate phase survives (see also the phase diagrams in 
Appendix [B]) . The dynamics of aggregates for large sys¬ 
tems size is studied in details in Sec. m 

The continuous transitions between aggregates and 
bands (or highly ordered phases) for N r-^ results from 
the competition between elastic energy and the impossi¬ 
bility of the system to sustain long-range polar order. For 
not too large system sizes, i.e. for N N^, the shape of 
the aggregates is roughly circular (Fig. 01 panel ii)) and 
at the center of the aggregate we find one single topo¬ 
logical defect: i.e. at the mesoscale, at the center of the 
aggregate we cannot define an average orientation for the 
rods. Due to the active forces, at the center of the aggre¬ 
gate rods are strongly compressed, which implies that the 
potentials Ui adopts high values (see Fig. 01 bottom row). 
This implies that when one of these aggregates is formed, 
the total elastic energy of the system Utot = Ylf=i Ui in¬ 
creases. On the contrary, in large polar structures such 
as bands, rods are roughly parallel to each other and 
therefore are much less compressed by their neighbors, 
and the total elastic energy is low. This is evident on the 
bottom left panel of Fig. 01 The dynamics at N 
can be summarized as follows. Large polar clusters form 
and eventually a band emerges, but since the system is 
too big for the band to remain stable, at some point the 
band breaks. The collapse of the band gives rise to the 
formation of new giant polar clusters which eventually 
collide head on leading to a large aggregate: a process 
reminiscent of a traffic jam. The formation of the aggre¬ 
gate leads to a sharp increase of the total elastic energy. 
Let us recall that forces and torques act in such a way 
that they tend to minimize Ui. In short, the system re¬ 
laxes by destroying the new formed aggregate, which give 
rise to the formation of new polar clusters and the cycle 
starts again. In larger system sizes, i.e. for N ^ N^, 
aggregates are more complex. This is addressed in the 
next section. 
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FIG. 5: Dynamics of an aggregate. Top panel: aggregate size 
and gas density as function of time. The aggregate bound¬ 
ary exhibits large fluctuations due to the emergence of ori¬ 
entational defects that lead to the detachment of large polar 
clusters from the aggregate. Second line: The three panels 
display in chronological order one of these events. The corre¬ 
sponding time window is indicated by the vertical grey area 
in the top left panel. Dashed circles indicate the detachment 
of polar clusters. The inset shows that topological defect. 
Bottom panel: Orientational and positional order within in 
the aggregates. Left: p{d) is the probability density that the 
distance between the centers of two rods is d (in units of i) . 
We observe peaks at multiples of the rod length i (the smaller 
peaks between 0 and 1 correspond to w). This is the finger¬ 
print of positional (smectic-like) order: the rods are arranged 
as on a lattice. Right: Snapshot of an aggregate. Inset: Zoom 
on a subdomain of the aggregate. The dashed black circle in 
the snapshot provides a visual estimate of what we call the 
core of the aggregate. Colors encode rod orientation to show 
the domains of similarly orientated rods. The arrows indicate 
the local orientation of the rods. Simulations of the aggre¬ 
gate dynamics correspond to A = 80000, rj = 0.15, and n — 7 
{N = 10000, g = 0.3, and = 4 for the aggregate in the 
bottom panel). For movies, see 0. 


V. AGGREGATE DYNAMICS 

Inside aggregates, the competition between active 
forces and local polar alignment leads to new physics un- 
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seen in other active systems. This is particularly evident 
for very large system sizes, i.e. for A/" ^ that is when 
aggregates are big enough to exhibit multiple topological 
defects of the local orientation of the rods. 

Let us recall that aggregates are formed by polar clus¬ 
ters of rods that are trapped inside these structures. The 
rods inside the aggregate do not only exhibit local orien¬ 
tational order, but also local positional order as indicated 
by the peaks exhibited by the probability density p{d) of 
finding the center of a rod at a distance d of the center 
of another rod (see bottom panel of Fig. [5]). Topological 
defects are areas where, at the mesoscale, as mentioned 
above, we cannot define an average orientation of the 
rods as illustrated in the bottom inset of Fig. [5] (areas 
where the arrows meet). In such areas, due to the active 
forces, rods are strongly compressed by the active push of 
all surrounding SPR (see inset in Fig. Since the com¬ 
pression is due to the presence of active forces, we refer 
to this phenomenon as active stresses. For N ^ we 
observe the emergence of multiple defects that lead to an 
increase of the elastic energy and the build-up of stresses. 
Notice that more topological defects imply larger values 
of the elastic energy. There are two clear consequence 
of the presence of multiple topological defects. On the 
one hand, aggregates are no longer roundish but rather 
irregular as illustrated in Fig. [5l On the other hand, 
now the system can relax the elastic energy by reduc¬ 
ing the number of topological defects. Notice that for 
N aggregates are relatively small and exhibit one 

topological defect, and thus, the only way to eliminate 
the topological defects is by destroying the aggregate. 
For A^ ^ A^*, given the presence of multiple topological 
defects, eliminating one topological defect does not re¬ 
quire to eliminate the aggregate. As a matter of fact, for 
very large system sizes, the interplay between topological 
defects and active stresses lead to large fluctuations of the 
aggregate boundary and aggregate mass (i.e. aggregate 
size) as indicated in Fig. [5l The most distinctive feature 
of the observed phenomenon is the large fluctuations ex¬ 
perienced by the aggregate mass correspond to ejections 
of remarkably large macroscopic polar clusters from the 
aggregate, that can be as large as 10% of the system size 
(i.e. involving more than 10^ rods). Fig. [5l By this pro¬ 
cess, i.e. the ejection of large polar clusters, the aggregate 
manages to decrease its elastic energy. The ejected po¬ 
lar clusters typically dissolve while moving through the 
gas phase outside the aggregate, leading to a sudden in¬ 
crease of the gas density, top panel in Fig.O This results 
in a higher absorption rate of SPR by the aggregate that 
starts again to increase its mass. The system dynamics 
- in the phase-separated state - can be summarized (in a 
simplified way) as follows: Aggregates grow in size and 
multiple defects emerge inside the aggregate. This results 
in active stresses that build up and lead to fluctuations of 
the aggregate boundary and ejection of huge polar clus¬ 
ters. This implies a reduction of the aggregate size and 
its elastic energy, and an increase of the gas density at 
which point the cycle starts again. 


VI. DISCUSSION 


The ejections of thousands of particles in densely 
packed and highly ordered clusters - the most distinc¬ 
tive feature of the described dynamical phase-separated 
state - requires the combined effect of an effective align¬ 
ment mechanism and active pushing (or stresses) acting 
among the particles. The combination of these two ele¬ 
ments is not present, to the best of our knowledge, in any 
other active matter model and is a distinctive property of 
SPR. And even for SPR such effects are only evident for 
large enough system sizes, i.e. above N^. For instance, 
self-propelled disks (SPD) exhibit active stresses but no 
alignment among the SPD While in SPD sys¬ 

tems phase separation is also observed, and high-density 
objects as aggregates are found, the dynamics of these ob¬ 
jects is totally different: SPD aggregates are not formed 
by polar clusters of SPD and they do not eject huge po¬ 
lar clusters of active particles. Moreover, SPD inside the 
aggregates are contained by a thin ring of SPD pointing 
inwards, which suggests that in the presence of several 
aggregates, particle exchange occurs via an evaporation¬ 
like process, controlled by the angular diffusion coeffi¬ 
cient, leading to a classical coarsening process 
which can be described by an effective Cahn-Hilliard 
equation 0. In sharp contrast with this scenario, the 
phase separation of SPR starts with a non-equilibrium 
ballistic clustering process that leads to the formation of 
polar clusters 0 , which in turn collide ballistically until 
eventually a traffic jam of polar clusters occurs and an 
aggregate emerges. The density instabilities and fluctu¬ 
ations observed with SPR are also remarkably different 
from what has been reported in systems of point-like ac¬ 
tive particles with an alignment mechanism as in the Vic- 
sek models HQ, where density fluctuations and phase 
separation (in the form of bands) require the presence of 
either polar or apolar long-range or quasi-long range or¬ 
der [^, which are absent in SPR as shown here. Finally, 
though in active particle systems that combine an align¬ 
ment mechanism and a density-dependent speed [10, [31| 
polar clusters, bands, and aggregates are found, as occurs 
in SPR for small system sizes, the large-scale properties 
of these systems are radically different from what has 
been reported here for SPR. In such systems, there is no 
active push among the active particles and while topo¬ 
logical defects cannot be ruled out, it can be safely stated 
that they cannot generate active stresses. Thus, in such 
models bands and polar phases might exist even in the 
thermodynamical limit, and certainly aggregates cannot 
eject polar clusters. 


In summary, the physics of SPR is remarkably different 
from the one observed in any other particle system. The 
coupling between active stresses and order - that implies 
that topological defects and active stresses are intimately 
related - leads to novel phenomena as the here reported 
phase-separated state characterized by the ejection of po¬ 
lar clusters. 
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Appendix A: Demonstration that Dp <C Da 

Here, we derive the diffusion coefficient for an isolated 
SPR subject to fluctuations. From the equations of mo¬ 
tion Eq. (1) and Eq. (2), it is clear that the temporal 
evolution of an isolated SPR can be expressed as: 



■\( 



gas 



clusters ^ 



unstable 



aggregate 




2 4 J5 8 10 0.075 0.15 „ 0.3 


= u(i) (Al) 

Oi = m, (A2) 

where in Eq. (EH) we have expressed explicitly in 
terms of the components parallel, V(6>^), and perpendicu¬ 
lar, to the long axis of the rod. The terms (t), 

and ^i{t) refer to independent delta-correlated 
noises [5l|, as in Eq. (1) and Eq. (2), but where we 
have absorbed the corresponding friction coefficients into 
the noise definition. The term u{t) in Eq. (jAip refers 
to the instantaneous velocity of the particle and should 
not be confused with the interaction potential. The po¬ 
sition of the SPR at time t can be formally expressed 
as = f^dsu(s), and the mean-square displacement, 
through the Taylor-Kubo formula [^, is given by: 

(xf(i))=/ ds f ds'{u{s).u{s')). (A3) 

Jo Jo 

To compute (u(s).u(s')) and finally to evaluate the inte¬ 
gral, we use that (d^||(t)) = {d-i±{t)) = = 0 and 

the autocorrelations of the independent delta-correlated 
noises [5l| that read: 

(<5'i||(s)<5-i||(s')) = 2L)||(5(s-s') (A4) 

(a-ij_(s)aj_L(s)) = 2Dj_S{s-s') (A5) 

(^i(s)li(s')) = 2DeS{s-s'), (A6) 

while other combinations vanish, i.e., (d^y = 

= 0. This correlations al¬ 
low us to compute (xf (t)), which for t ^ takes the 

form: 

(x2(t))=2At + 2(D||+Dx)i. (A7) 

D'e 

The diffusion coefficient is defined in two dimensions as 
D = limt^oo(xf (t))/(4t), and thus we obtain: 


FIG. 6: Exploration of the parameter space (phase diagram). 
Each squared point in the plots correspond to a parame¬ 
ter set {A, ? 7 , Ac} for which we have performed simulations. 
The color of a point encodes the observed state of active 
matter. Black: disordered, not phase-separated state (gas). 
Red: ordered phase-separated state (polar clusters or polar 
band). Green: disordered phase-separated state (aggregate). 
Blue: unstable situation, where the system oscillates between 
the ordered phase-separated state and the disordered phase- 
separated state. The color code for the regions of the phase 
diagram is the same as for the simulation points. 


the ratio between the corresponding drag coefficients: as 
example D\\/D± = C-l/C||- From this, we learn that while 
Dp (X T, Da oc 1/T. Einally, we notice that we can ex¬ 
press, for the reasons given above, Dy = Dq(o/(\\ and 

D± = Do(q/(±. Eor the values of Dq^ uq, C ||5 C-L Ce 
used in the main text, the ratio between the active and 
passive diffusion coefficient is such that DpjDa ^ 10“^, 
i.e. Dp Da. In such a regime, we can ignore the contri¬ 
bution of Dp, meaning that we can neglect d^y and 
and focus exclusively on the role of on Da- 

Appendix B: Exploration of the parameter space 

The left-hand panel of Eig. [6] indicates, the parameter 
sets {A^, /i:} for which we have performed simulations at 
packing fraction = 0.3, and the observed states of active 
matter. Note that when increasing N, the orientationally 
ordered phase {i.e. clusters) progressively disappears and 
is replaced by an disordered phase (aggregate). 

The right-hand panel of Eig. [6] corresponds to the phase 
diagram for large systems {N = 80000). We observe that 
the area of the phase space corresponding to orientational 
order (clusters) is reduced to a small band (which should 
completely disappear when further increasing N). 

Appendix C: Polar order and clusters 


D = 


2De 


+ 


D\\ + D± 


— Da + D 


P 5 


(A8) 


where we define the active diffusion coefficient by Da = 
and the passive diffusion coefficient by Dp = ^ 

Notice that if the noise terms are of thermal origin 
in Eqs. EH and (jA2p . then D\\ (x T, D±^ (x T, and 
Dq oc T, where T is the temperature. The ratio between 
any of these coefficient, e.g. D\\/D±_ is proportional to 


In Eigs. 1 and 3 of the paper we have investigated 
the global polar order Si in the system. We observed 
polar order for /i: > in Eig. 1 and for A < in 
Eig. 3. If there are only few giant polar clusters (or 
even a single percolating band) in the system, the ob¬ 
served global polar order directly results from the polar 
order within these few polar clusters. However, if there 
are many clusters, their orientations may be correlated, 
which also contributes to global polar order. Eor very 
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FIG. 7: Left: cluster size distributions p{m) and the 

corresponding polar-order-within-clusters-distributions *S'i,m 
(i.e. the polar order within clusters as function of their 
size m). Right: Global polar order parameter effectively ob¬ 
served Si and global polar order parameter obtained under 
the hypothesis that cluster orientations are fully decorrelated 
Si,uncorrelated- Simulations correspond to 77 = 0.15, = 10 

(and N = 20000 for the iSi^m-distribution, but this distribu¬ 
tion is very similar for the other values of N). 


big systems {i.e. in the thermodynamic limit), global 
polar order can only be observed if the orientations of 
the clusters are correlated. Here we attempt a quantifi¬ 
cation of the impact of cluster correlations on global po¬ 
lar order. Therefore, we assume that polar order within 
clusters only depends on cluster size m and denote it 
as Si^m = ^lnmJ2h=i the polar order 

associated to the cluster of mass m. We then measure 
the cluster size distribution p{m). This permits us to 
compute the value of the global order polar parameter 
Si,uncorrelated uudcr the hypothcsis that there are no cor¬ 
relations, i.e. we assume that, at a given time t, any 
cluster j has a moving direction Oj {t) (defined as the av¬ 
erage moving direction of all the rods within the cluster) 
which is distributed according to a uniform random dis¬ 
tribution. 

The mathematical derivation is the following: Starting 
from the definition of the global polar order parameter 
given by Eq. 3 of the main text, we rewrite the average 
over all particles as an average over all clusters, i.e. 


where rrij is taken from the steady state cluster size 
distribution p{m) and Oj from an homogeneous distribu¬ 
tion between [0, 27r], such that the number of rods is N. 
This is an approximation because we replace the polar 
order of a cluster by the average order of a cluster of the 
corresponding mass, The distributions p{m) and 

Si^rn have been directly measured from simulations, see 
left and central plot of Fig. [71 The final step is to assume 
that each cluster points in a (uniform) random direc¬ 
tion 0j{t) as explained above. The computation of the 
global order parameter is done by performing a Monte 
Carlo algorithm, required to average over the random 
direction of the cluster (and the random cluster size). 
We refer to this quantity as Si^uneorrelated, see Fig. [71 
We then compare Si^uneorrelated to the effectively ob¬ 
served polar order Si (see the right-hand plot in Fig. [T]). 
For all investigated packing fraction p (but in particular 
for high packing fractions) and system sizes N, the ob¬ 
served polar order (Si) is systematically bigger than the 
one resulting from the cluster decorrelation hypothesis 
{Si,uneorrelated)- This indicates the presence of cluster 
correlations. Nevertheless, we observe that Si becomes 
closer to Si^uneorrelated whcu the systcm size N increases. 
This suggests that cluster correlations are of finite size 
(i.e. the orientations of clusters which are sufficiently far 
away from each other are no more correlated). This is 
a strong indication that in the thermodynamical, limit 
cluster-cluster correlations are weak and cannot produce 
(global) polar order. 


^1, 


uneorrelated 



N 
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